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Continuous time random walks are non-Marlcovian stochastic processes, which are only partly 
characterized by single-time probability distributions. We derive a closed evolution equation for 
joint two-point probability density functions of a subdiffusive continuous time random walk, which 
can be considered as a generalization of the known single-time fractional diffusion equation to two- 
time probability distributions. The solution of this generalized diffusion equation is given as an 
integral transformation of the probability distribution of an ordinary diffusion process, where the 
integral kernel is generated by an inverse Levy stable process. Explicit expressions for the two time 
moments of a diffusion process are given, which could be readily compared with the ones determined 
from experiments. 
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The concept of continuous time random walks (CTRWs), introduced by MontroU and Weiss [IJ almost four decades 
ago, has been successfully applied to a wide variety of transport problems in physics 2 . In recent years the relationship 
between CTRWs and a class of Fokker-Planck equations with fractional temporal and/or spatial derivatives has 
attracted a lot of attention [5]. In this line of research a major focus has been on the investigation of anomalous 
diffusion T, a nonequilibrium phenomenon occurring in scientific fields ranging from astrophysics to biophysics and 
econophysics. A well known example is a subdiffusive CTRW in a force field, which leads to the well-known fractional 
Fokker-Planck equation (FFPE) [5 , invoking the fractional Riemann-Liouville differential operator [6j. This equation 
can be written in the equivalent form: 

oC^f{x,t) = LFpf{x,t), (1) 
where qC^ denotes the Caputo fractional differential operator [B], defined as 

oC^git) ^ — ^ f g^'\t'){t-tr" (2) 
r(i - a) Jo 

for < Q! < 1 (g'^^^t) is the first derivative of g(t)). Lpp is a Fokker-Planck operator. For simple diffusion processes 
in Id Lpp — Dd^/dx^ and Eq. ([l]) is referred to as fractional diffusion equation (FDE) O |7]. In the limit a — > 1 
Eq. ([l]) reduces to the ordinary Fokker-Planck equation. Clearly, the Caputo fractional time derivative expresses the 
non-Markovian character of the underlying CTRW: the probability of finding the random walker at point x at time t 
depends on the whole history of the process from time up to t. In combination with fractional evolution equations, 
the CTRW constitutes a versatile stochastic model which can take into account sub- and superdiffusive behaviour 
in complex systems as diverse as turbulence [H], optical lattices, and biological cell motility (see [S] for an extensive 
review). Recent new developments concern a more fundamental understanding of CTRWs and their application, e.g. 
the connection to ageing phenomena [TOj , ergodicity breaking |IT] and inertial particles diffusing in a potential jl2j . 

However one important aspect has long been neglected in the literature on CTRWs and FFPEs. In general the 
ordinary Fokker-Planck equation is a deterministic evolution equation for the transition probability of the stochastic 
process, which is assumed to be Markovian. By virtue of the Markovian property one can calculate arbitrary ri-point 
pdfs with a single transition probability. On the other hand, in the case of the fractional analog Eq. ([T|, single-time 
pdfs are determined, which contain only very limited information about the stochastic process. In order to obtain 
complete information about a non-Markovian process like the CTRW, one has to consider the infinite set of n-point 
pdfs. Recently, such multi-point statistics have been investigated in different approaches |13l I14j . In |13j . a two-point 
Green function has been determined with the help of the FFPE backward propagator and allows for the computation of 
three-point fluorescence lifetime correlation functions. Our previous work [T3] discussed joint probability distributions 
of a CTRW based on a representation in terms of coupled Langevin equations. The purpose of the present paper is to 
derive a closed evolution equation for the two-point pdf of a CTRW, which is independent of the inverse Levy-stable 
process investigated in jT4]. This equation is the generalization of the single-time FDE ([T]) to two times. Considering 
the widespread interest in Eq. ([l]), the two-point generalization is expected to provide important further insight into 
the application of CTRWs to real world systems. In fact, multi-point statistics have already been probed in recent 
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experiments on protein conformational dynamics |15l 116] and highlight the importance of a more complete theoretical 
understanding of these statistics in anomalous diffusive systems. 

A suitable starting point for the investigation of CTRWs is a representation in terms of coupled Langevin equations 
introduced by Fogedby \17\. Here the motion of a Brownian particle in an external force field F{x) in d = I dimensions 
is described as: 



dX{s) 
ds 

dt{s) 
ds 



F{X) + r^{s), (3) 
r{s). (4) 



In this framework the CTRW is parametrized by the continuous path variable s, which may be regarded as arclength 
along the trajectory. X denotes the physical space and t the physical time. Both are given as stochastic processes 
in the 'eigentime' s. Their statistics are determined by the properties of the stochastic variables 77(5) and t(s). In 
this work we consider the special case of statistically independent increments rj and r, i.e. jump lengths and waiting 
times are uncoupled. Furthermore we want to restrict our considerations to a subdiffusive CTRW without force field. 
Accordingly, ry(s) is assumed as a standard Langevin force with properties (?/(s)) = and {ri{s)r](s')) = C^is — s') . 
X(s) is consequently given as a standard Markovian Wiener process. The subdiffusive characteristics enter via the 
process t{s). Its increments r(s) are assumed to be broadly distributed such that t{s) constitutes an asymmetric 
Levy-stable process of order a with < a < 1. Levy-stable processes of this kind induce a diverging characteristic 
waiting time {t{s)). We are interested in the process X{t) — X{s{t)), i.e. the behaviour of the physical space variable 
X as a function of physical time t. The pdfs of this process are defined as (e.g. for two points): 

f{x2. t2; Xi,h) = {Six2 ~ X(S2))<5(S2 - s{t2))Kx^ " - s(ti))) . (5) 

As a consequence of above specifications, X(i) is non-Markovian and reveals subdiffusive characteristics: (a;(t)^) ~ t". 
Its statistical properties are closely related to the properties of the process s(i), the inverse of the process i(s). The 
investigation of the statistics of the inverse Levy-stable process s(i) has been the focus of reference [M]. One of the 
main results is the following integral transformation for the n-point pdf /({x^jti}) of the process X{f) [T3]: 

!{.{xiMY) = / dsi... dsn h{{si,ti})fM{{xt,Si}). (6) 
Jo Ja 

h{{si,ti}) denotes the n-point pdf of the process s(t) and fM{{xi, Si}) the n-point pdf of the Markovian process X(s). 
Eq. ([g]) states that the pdf of the non-Markovian process can be determined by a transformation of the corresponding 
Markovian process. The integral kernel is generated by the inverse Levy-stable process s{t). Whereas the pdf 
fMi{xi, Si}) is obtained in a straightforward way from the Langevin equation ([3|, the determination of h{{si,ti}) 
proves to be the crucial point. In the one- and two-point case h assumes a simple form in Laplace space. 

In the following we present the derivation of the two-point fractional diffusion equation starting from the trans- 
formation Eq. ([6|. The derivation of the single-time FDE (|T| proceeds along the same lines and is not shown here. 
We define Laplace transforms in the usual way: C{g{t)} = g{t)e~^*dt, accordingly for multiple variables ti A^. 
Functions with argument A^ always denote the Laplace transform if not otherwise indicated: g{X) = C{g(t)}. In 
Laplace space Eq. Q for the two-point pdf f{x2,t2',xi,ti) assumes the similar expression: 



/(x2, A2;a;i, Ai) = / dsi ds2 h{s2,\2\si,\i)jM{x2,S2\xi,si). (7) 
Jo Jo 

The Laplace transform of the two-point pdf h is explicitly given as [TIJ : 

his2,X2;s„X,) = ^(.2-.i) ^°"^^\^^^^° + ^^ e-^(^^+^-)° 

MM 

+e{s2 - ,^)Mi(^i±M!l^Me-(Ai+A.r«ie-A?(s.-s,) 
Ai A2 

+e(si - S2) t ~ ^?) e-(A.+A.r..^-A?(.,-s.)_ 

A1A2 

An equation for h{s2, A2; si, Ai) is readily derived in the form 

d d X 

7^ + 7^ ) h{s2, A2; Si, Ai) = -(Ai + A2)"/i(s2, A2; Si, Ai), (9) 
dsi ds2 / 
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with initial conditions 



M0,A2;0,Ai) - 



A? - (Ai + A2)° + 
A1A2 



MM 



MO, A; 



,.„A0 = (A?)((A. + A2r-A?)^_,,, 



A1A2 



(10) 



In general the two-point pdf /M(a;2) ^2; Xi, Si) of the Markovian diffusion process is determined by the transition 
probability: fM{x2,S2',xi,si) = P{x2, S2\xi, si)f{xi, si). In turn P(a;2, S2|a;i, si) is obtained from the ordinary 
Fokker-Planck equation. Due to the Markovian property it is not necessary to formulate a closed evolution equation for 
fM{x2, S2',xi, Si), however such an equation can be derived in a straightforward way from the two-point characteristic 
function of the process X{s). The result is: 



_d_ _d_ 
dsi ds2 



fM{x2,S2;Xi,Si) = L{X2, Xi) fM{x2, S2; Xi, Si), 



where we define a generalized diffusion operator as the inverse Fourier transform 

L{x2,xi)g{x2,xi) 



(11) 



(12) 



2 \dxl dxidx2 



dxl 



g{x2,xi). 



The definition Eq. ( |12| holds in the same form in higher dimensions. Furthermore, broadly distributed symmetric 
increments rj{s) in Eq. (|3]) would result in a characteristic exponent < /Lt < 2. The inverse Fourier transform then 
leads to a two-point generalization of the fractional Riesz/Weyl operator The initial conditions of Eq. (Ill read: 



fM{x2,S2\Xi,Si = {)) = fM{x2,S2)S{xi), 

.fM{x2,s2 = 0;xi,si) = /M(a;i,si)i5(x2), 
/M(a;2,S2 = 0;a;i,si = 0) = S{x2)S{xi). 



(13) 



Now we can proceed as follows. We multiply Eq. ^ by (Ai + A2)" and substitute Eq. After performing partial 
integrations with respect to si and S2, Eq. (11) can be substituted leading to the following equation: 



(14) 



The specific boundary terms occur due to the partial integrations and assuming h{s2 ~^ 00, t2',Si, ti) — h{s2, ^2! Si ^ 
00, ti) = . These terms can be absorbed into the Caputo fractional differential operator, which we generalize to two 
times in a straightforward way. First we note that h{s2, A2; si — 0, Ai) and h{s2 = 0, A2; si, Ai) are determined from 



(Ai-|-A2)"/(x2,A2;xi,Ai) = i(a::2, xi)/(a;2, A2; a;i, Ai) + / ds2 /i(s2, A2; 0, Ai)/Af(a;2, S2)5(a;i) 

JO 

poo 

+ / dsi /i(0,A2;si, Ai)/A/(xi,si)(5(a::2). 



Eq. (If: 



/i(s2,A2;0, Ai) = S{s2) 



A? - (Ai + A2)" + A£ ^ ^^^^ _^ _ 1 ^.-i^-A=.. 

Ai 



A1A2 



(15) 



h{s2 — 0, A2j_si, Ai) takes the analogue form. Both expressions can be substituted into Eq. (14 1. In the second 
term of Eq. (15) the Laplace transform of the single-time distribution h{s2,t2) occurs [H]: h{s2,X2) = A" "^e"'^?''" . 
Consequently, with the help of Eq. ^ for the case n = 1 the integrations in the boundary terms of Eq. ( 14 1 lead to 
the single-time distributions /(a:2,A2) and /(a;i,Ai): 

(Ai + A2)"/(a;2, A2; Xi, Ai) = L{x2,xi)f{x2,\2\ xi,\i) + ((Ai + A2)" - A^)^/(a;2, \2)5{xi) 

M 

+((Ai + A2)" - X^i)^f{xi,Xi)5{x2) + 2 ^°~^^\+^^"^° + ^^ ^(x2)J(xi). (16) 
M MM 
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Rearranging terms yields: 

(Ai + A2)^ 



A1A2 



[AiA2/(a;2, A2; xi, Ai) - X2)S{xi) - Xif{xi, Xi)S{x2) + S{x2)S{xi)] 



+ [^^f{x2, A2) - \r'S{x2)f-^ + [A?/(xi, Ai) - ArM(a:i)]^ 

Ai A2 

L{x2,Xi)f{x2,X2;Xi,Xi) + — (^l + -^2) + ^2 ^^^^-jj^^^-^^ 

Al A2 



(17) 



In order to perform the inverse Laplace transform of this equation we state the two results (r(/3) denotes the Gamma 
function) : 



c^eit2-h) 



r(i-a) 



Q{h-t2) 



r(i-a) 



(Ai + A2)^ 
A1A2 



and 



C 



d_d_ 
dhdfi 



g{t2, h) )■ = A2Ai<?(A2, Ai) - A2g(A2, ti = 0) - Ai,g(t2 = 0, Ai) + g{t2 ^0,h = 0). 



(18) 



(19) 



Furthermore, due to the definition Eq. ([5]) the boundary conditions of the pdf f{x2, ^2; xi, ti) are as in the Markovian 
case Eq. (13 1: f{x2,t2;xi,ti = 0) = f{X2,t2)6{xi), etc. The Laplace transform of the single-time Caputo operator 
Eq. ^ reads C{oC"g(t)} = X'^g^X) — X"~^g(t = 0). Having all this in mind, the inverse Laplace transform of Eq. (17 1 
then yields our main result: 



0C^^,J{x2,t2;Xi,h) + oC^J{x2,t2)S{Xi) + oC^J{xut,)d{x2) 



Lix2,Xi)f{x2,t2; Xi,ti) 

+h{0,t2;0,ti)S{x2)6{xi). 



(20) 



This is the two-point fractional diffusion equation, generalizing the single-time FDE to two-point probability distri- 
butions. Here, a generalization of the single-time Caputo operator to two times is introduced as: 



oCr,,,g(i2,ii) 



1 



- [6(^2 -ii)<r" + 0(^1 -*2)i^"] 



d d , 



(21) 



r(l-a) 

where the ** denote a double Laplace convolution with respect to ^2 and ti. With the definition of the two-time 
fractional Caputo derivative Eq. (21 1 the general form of the diffusion equation is clearly visible. In addition, 
boundary terms occur which can be interpreted according to the underlying random walk. QC"^f{x2,t2)S{xi) and 
oC"^f{xi,ti)S{x2) describe the propagation when either ti = or t2 = 0. h{0,t2;0,ti)6{x2)S{xi) is due to the non- 
zero probability that the random walker stays at the initial site. The analytical form of this term is obtained by 
performing the inverse Laplace transform of h{0, A2; 0, Ai) (Eqs. ( 10 1): 



Mo,i2;0,ti) = e{t2-h)Ho,t2) + e{h^t2)Ho,h), 



(22) 



where h(0, ti) = t^" /T(l — a). This expression agrees with the interpretation given above. For t2 > ti, the probability 
of staying at the initial site until t2 is given by the single-time distribution h{0, ^2) and vice versa. As in the 
single-time case, the two-point FDE reduces to the form of its Markovian counterpart Eq. (Ill for a ^ 1 (see the 
expression in Laplace space Eq. (17l). Also, the transformation Eq. ^ reproduces this result, since for a ^ 1: 
h{s2,t2', si,ti) = S{s2 — t2)S{si — t^{see Eq. (|8|). The occurrence of additional boundary terms in Eq. (20 1 can be 
considered as a general signature of non-Markovian processes and is even more prominent in the multiple-time case. 
A generalization to n-point probability distributions is obtained in a straightforward way along the lines outlined 
above. 

From Eq. ( 20 1 the two-time moments of the subdiffusive CTRW can be calculated without invoking the properties 
of the inverse Levy-stable process s{t). Consider e.g. the simplest moment {x{t2)x{ti)) . Multiplicating Eq. (20l with 
X2, xi and integrating from —00 to 00 yields oC't"ti {x{t2)x{ti)) — where C is the strength of the stochastic force 
in the Langevin equation ([3]). A solution of this equation can be calculated in Laplace space. Since we use the usual 
convention X{t = 0) = the boundary terms due to the Laplace transformation of oC^^ti all vanish. One immediately 
obtains the solution C{{x{t2)x{ti))} = C/(AiA2(Ai -t- A2)"), which reads in real space: 



{x{t2)x{h)) = Q{t2~h)C 



4-a 



+ eih-t2)C 



r(a-t-i)' 



(23) 



Time tj Time 

FIG. 1: Semi-logarithmic plot of the correlations oi x{ti), x{t2) and x{ti)^, x{t2)'^ ■ Here, ti is fixed at 100, setting the symmetry 
axis, and C — 1- units are non-dimensional. 



Due to the fact, that f{x2,t2',xi,ti) is given as the transformation Eq. (|6| only even higher order moments are 
non-zero. They obey the general recursion formula (for non- negative integers m, n): 

oCZu {xihf"'xihf") - m{2m - 1)C (x(t2)"""'x(ti)2") + n{2n ~ !)( {xihf'^xihf"-^) 

+4TOnC (x(i2)"""'a;(ii)2"-i) . (24) 



In the case m = n = 1 we obtain: 



-a; a 



t2 



Here, F{a,b; c; z) denotes the hypergeometric function (see e.g. [IH])- Figure [l] shows a semi-logarithmic plot of the 
correlations corresponding to Eq. (23 1 and Eq. (251 as functions of ^2 for three different a- values. Both correlations 



exhibit a clear power law decay for ^2 S> ii: 0:^2 cx i^" respectively. It would be interesting to determine 

these correlations from experiments. 

We have demonstrated that a consistent generalization of the well-known single-time fractional diffusion equation to 
two-point pdfs can be derived on the basis of the coupled Langevin equations introduced by Fogedby as a representation 
of CTRWs. Special features of the two-point FDE are a two-time fractional differential operator of the Caputo-type 
and the occurrence of additional boundary terms. Its solution is expressed in terms of an integral transformation of 
the two-point pdf of the corresponding normal diffusion process. As in the single-time case the a ^ 1 limit reduces 
the FDE to its Markovian counterpart. Furthermore we derived recursion relations for arbitrary two-time moments 
of the subdiffusive CTRW. It should be noted that the occurrence of fractional derivatives is a consequence of the 
properties of the stochastic process t(s) which determines the temporal behaviour of the CTRW. Thus the derivation 
of evolution equations, as presented in this paper, should apply to a whole class of systems, which can be described 
by two independent stochastic processes for X and t. Here the simplest case has been solved, namely X specified as 
a Wiener process. For this case, we have determined explicit expressions for the two-time moments, which could be 
readily compared with moments obtained from experiments. An investigation of the multi-point statistics of other 
CTRW related processes, e.g. the anomalous diffusion of weakly damped inertial particles [H], is left for future work. 
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